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COMPUTER PROGRAM FOR CALCULATING ISOTHERMAL, 

TURBULENT JET MIXING OF TWO GASES 
by Leo F. Donovan and Carroll A. Todd 
Lewis Research Center 

SUMMARY 

Fluid mechanics is perhaps the most significant area in the investigation of the 
coaxial-flow, gas-core nuclear reactor. This report describes a computer program for 
solving a simplified model of the turbulent mixing that occurs between a central fuel jet 
and a surrounding, faster-moving coaxial stream of propellant. As such, this report 
constitutes a step toward a better understanding of this aspect of the gas -core nuclear 
rocket. Local values of time -averaged velocity and the mass fraction of fuel can be 
calculated for a reactor in less than a minute on the IBM 7094. 

The von Mises transformation was used to convert the axisymmetric forms of the 
isothermal boundary layer momentum and diffusion equations to forms amenable to nu- 
merical solution. The effects of confining walls were not considered. The program can 
solve problems in which the initial velocities and densities of the two streams differ 
greatly, by using expressions for eddy viscosity that vary radially as well as axially. 

The effects of initial coaxial-stream- to jet- velocity and density ratios on velocity 
and mass fraction profiles are shown. An unspecified reference density occurs in the 
eddy viscosity formulation; the influences on velocity and mass fraction profiles of 
several choices for the reference density are illustrated. Radial and axial variations of 
eddy viscosity and the product of density and eddy viscosity are also given. Estimates 
are made of the effects of initial velocity ratio and initial density ratio on the mass of the 
major jet component contained within a given volume. The maximum and minimum 
amounts that could be present are included for comparison. 


INTRODUCTION 

The concept of a coaxial-flow, gaseous nuclear reactor provides the motivation for 
the work to be described in this report. A brief discussion of some of the features of 
this concept will be helpful in understanding the relevance of turbulent jet mixing. 



An important characteristic of rocket performance is specific impulse; specific 
impulse is the thrust obtained per unit weight flow rate of propellant mixture expelled 
from the engine. High specific impulse is desirable since less propellant will be required 
to produce a given total impulse (i. e. , integral of thrust over the operating time of the 
reactor) . Specific impulse is approximately proportional to the square root of the ratio 
of exhaust temperature to molecular weight; thus , high temperatures and low molecular 
weights are desirable. Chemical rocket performance is limited by the temperature to 
which the heat of combustion will raise a given fuel-oxidizer combination; for example, 
advanced hydrogen-oxygen rockets produce a specific impulse of about 450 seconds. The 
great advantage of nuclear rockets is the high specific impulse that can be obtained by 
using hydrogen as the propellant. Solid-core nuclear reactors, however, must operate 
at temperatures that fuel-bearing materials can withstand and are thus limited to specific 
impulses of about 1000 seconds. 

Higher fuel temperatures are possible in the gas-core nuclear reactor since the fuel 
is not supported on solid surfaces. Rather, a slowly moving gaseous fuel mass radiates 
heat to a coflowing annular propellant stream. The initial ratio of propellant to fuel 
velocity must be high in order to keep the loss of fuel as low as possible. With such a 
reactor, specific impulses of 2000 to 3000 seconds may be possible. Use of nuclear fuel 
and hydrogen propellant results in small initial propellant- to fuel-density ratios. 

The gas -core nuclear reactor problem was made more amenable to solution by 
dividing it into three major areas (ref. 1) (viz, nuclear aspects, radiant heat transfer, 
and fluid mechanics). The goal is to recombine the separate parts into a meaningful 
whole after each relevant process is understood. The first two areas have been discussed 
in part elsewhere (refs. 2 and 3); fluid mechanics, including mass transfer, was con- 
sidered in reference 4 and is treated in the present report. In reference 4 molecular 
transport coefficients were retained, and eddy viscosity was assumed to be a function of 
axial position raised to an arbitrary power which was determined by comparison with 
data from a bromine -air experiment. 

In the present analysis a modified Prandtl eddy viscosity was used along with the fact 
that in turbulent jet mixing the molecular transport coefficients are negligible compared 
with the turbulent transport coefficients. The equations governing the mixing are the con- 
tinuity equation, the momentum equation, and the diffusion equation. Boundary layer as- 
sumptions were used to simplify these equations. A sketch of the model analyzed is shown 
in figure 1. The problem at hand differs from what has been solved before in that a large 
initial coaxial -stream- to jet-velocity ratio is coupled with a small coaxial-stream- to 
jet-density ratio, so that an appropriate formulation of the eddy viscosity is not known. 

The problem was simplified for this report by eliminating the effect of confining 
walls and considering a free jet. Most of the work in jet mixing has been done on free 
jets, and the results (ref. 5) of these experiments and analyses can be used as a basis 
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for estimating eddy viscosity. Prandtl's postulate that eddy viscosity is proportional to 
the product of half radius and maximum velocity difference has been shown to agree with 
data in the far jet region, when density is constant and the initial coaxial-stream- to 
jet-velocity ratio is small. However, there is little agreement on how this postulate 
should be modified when density varies. Ting and Libby (ref. 6) have extended this 
formulation to allow for density differences, but adequate experimental verification has 
not been obtained. 

This report describes a computer program that has been developed for isothermal, 
turbulent jet mixing of two gases. The program is capable of solving problems with 
large initial velocity and density ratios without the restrictive assumption of constant 
eddy viscosity. The effects of velocity and density ratios on the mean flow properties 
are shown; also, the marked influences of the reference density, unspecified in the 
Ting -Libby formulation, are illustrated. 


ANALYSIS 

Boundary Layer Equations, Initial and Boundary Conditions 


The equations that are used to describe the isothermal, turbulent jet mixing of two 
gases are the time -averaged continuity equation and boundary layer forms of the momen- 
tum and diffusion equations. The von Mises transformation converts the momentum and 
diffusion equations to forms that satisfy the continuity equation identically. Eddy 
viscosity is specified empirically by a combination of Prandtl's constant density formu- 
lation and a relation proposed by Ting and Libby (ref. 6). 

For jet mixing at large Reynolds number, molecular transport is negligible com- 
pared with turbulent transport and can be ignored. At low Mach number, density changes 
result solely from mixing, and for a free jet the pressure is constant. Using capital 
letters to denote dimensional quantities, the axisymmetric forms of the continuity, 
momentum, and diffusion equations are as follows (ref. 5): 


— (PUR) + — (PVR) = 0 (1) 

ax aR 


PU — + PV — = — — (PER ' 
ax 3R R 3R V 3R ; 


( 2 ) 
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(Symbols are defined in appendix A.) The diffusion equation is written in terms of the 
mass fraction of component 1, the major (or total) constituent of the initial jet. It is 
assumed that the turbulent momentum and mass fluxes can be represented as the product 
of an eddy viscosity and the gradient of a time -averaged quantity. In addition, the 
turbulent Schmidt number Sc t is used to relate the eddy viscosities for momentum 
and mass. 

The initial and boundary conditions for the problem are as follows : 


u = u r 
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U — U e , Y — Y g R - co X 0 (6) 

These conditions may have to be modified when comparing computer results and experi- 
mental data. If the wall thickness of the jet discharge tube is not small compared with 
the tube diameter, the presence of a wall may significantly influence the early develop- 
ment of the flow. Also, the jet and coaxial-stream velocities will not, in general, be 
uniform but will have some distribution. If the duct surrounding the coaxial stream is 
not large compared with the jet discharge tube, the assumption of a coaxial stream of 
infinite extent is not justified. 

The equations can be made dimensionless in terms of the initial jet velocity, mass 
fraction, density, and radius. When lower-case letters are used to denote dimensionless 
quantities, the equations become 
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where 


pe = 


PE 


Wi 


The dimensionless initial and boundary conditions are as follows: 


u = 1, y = 1 0 < r < 1 

u = u e , y = y e r > 1 


X = 
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lH=0, v= 0,^=0 r = 0 x > 0 (12) 
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The density and mass fraction of component 1 can be related by the ideal gas law. 
Thus, 

m. 


P = 


M 


i 0 - 1 + — 

2 Yi 


M j (m 9 - 1) y + — 


(14) 


where mg = Mg/Mj. For most applications, the initial jet will be pure component 1, 
and the initial coaxial stream will be pure component 2. Then, p g = mg and the mole 
fraction of component 1 is 


c = 


1 - Pc 


(15) 


A "compatibility condition" that must be satisfied by the numerical solution can be 
obtained by evaluating the momentum equation at the centerline; this condition provides 
a check on mesh size. Thus, when 1’ Hospital's rule is used to evaluate the indeterminate 
form that arises, 
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von Mises Transformation 

The numerical integration can be facilitated by using the von Mises transformation 
(ref. 7, p. 136) to convert from (r,x) to (ip, x) coordinates, in which the continuity equa- 
tion is satisfied identically. Defining a stream function ip such that 
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and using the chain rule for differentiation 
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results in the following forms of the momentum and diffusion equations: 
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In these coordinates the initial and boundary conditions are as follows: 
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U = 1, y = 1 
U = u e , y = y e 

= o, = 0 

d\p 

u = V y = y e 

The transformation obtained by integrating 
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he first of equations (17) 






pur dr 


(25) 


can be used to determine \f/. 


Eddy Viscosity 


It remains to specify the eddy viscosity. Two formulations are required: one for 
the "near" jet, before the centerline velocity begins to change, where the mixing is 
more nearly planar; and the other for the "far" jet, where the mixing is truly axisym- 
metric. No attempt was made to make the eddy viscosity continuous at the point where 
one formulation replaces the other. 

Ting and Libby (ref. 6) have postulated relations between the eddy viscosities in 
constant -density and variable -density flows. These relations can be written as 


e = 



£* 


(26) 


for the near jet and 

e = e* f—\ — I 2 -£- r dr (27) 

\ p ) r 2 «/q p o 

for the far jet. The asterisk refers to constant-density flows and p Q is a reference 
density for the flow. Since the reference density is not specified, it must be determined 
by comparison of calculation and experiment. The value for the centerline eddy viscosity 
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can be obtained by expanding p in a Taylor series in r about r = 0, performing the 
integration, and taking the limit as r — 0. In this way 


€ t = 



( 28 ) 


The constant-density eddy viscosity can be represented by the following equations 
(ref. 7) in the near and far jets, respectively, 

e* = kjxfag - 1) (29) 

and 

e * = k 2 r l/2 (u e “ (30) 


For a jet discharging into a quiescent ambient stream, kg was found experimentally to 
be 0. 0256 (ref. 7). In terms of these expressions for eddy viscosity, the centerline 
compatibility conditions for the near and far jets, respectively, become 


and 



(31) 




du. 


dx 


= 2k 2 r 1/2 ^ u e ‘ 




(32) 


The unspecified reference density is presumably the centerline density, the coaxial- 
stream density, or a combination of these. It was assumed that a simple linear combina- 
tion was adequate. The reference density was thus taken to be 


p 0 = A p t + Bp e (33) 

where A and B are positive input constants. Restricting the sum of A and B to 1 
bounds p Q between p^ and p 0 . 
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Calculations 


In addition to the velocity, mass fraction, and mole fraction variations, several 
other quantities are of interest. The turbulent momentum and mass fluxes are 


3u 

- pe- 

(34) 

ar 


_ p£_ ii 

(35) 

SCj. 3r 



The width of the jet is characterized by its half radius (i. e. , the point at which the local 
velocity equals the average of the centerline and coaxial-stream velocities). 

A gross measure of concentration that can be easily obtained experimentally is the 
amount of light attenuated by an opaque substance. The attenuation can be related to 
concentration if the optical properties of the components are known. This technique was 
used in reference 8 with a bromine jet and a coaxial stream of air. This measure of 
concentration can be calculated as follows: 


c* 



P-Pe 
1 -Pe 


dr 


(36) 


Gas-Core Nuclear Reactor Calculations 

For these calculations, it was assumed that the fuel instantly vaporizes upon enter- 
ing the reactor. Heat-transfer calculations (ref. 3) were used to estimate average fuel 
and propellant temperatures so that molecular weights and densities could be calculated. 

The major fluid-mechanical figure of merit for the gas -core nuclear reactor is 
the amount of fuel contained within a given volume 

W = 2ttP.Y.R?I (37) 

where 

rl f d/2 

I = / / pyr dr dx (38) 

•A) * / 0 
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The calculations were performed for a coaxial stream infinite in extent , whereas in a 
gaseous nuclear rocket the coaxial stream is, of course, bounded. Indeed, it is now 
thought that the reactor radius will be only about twice the radius of the jet discharge 
tube. However, the radial integration in equation (38) was terminated at the position 
corresponding to a reactor wall in order to estimate the fuel content. Since the radial 
velocity at this position is not zero in the calculations, the situation to which the calcu- 
lations apply corresponds to a porous -wall reactor with this distribution of radial 
velocity at the wall. Alternatively, if the integration is carried out to a constant \p 
value corresponding to the initial reactor radius, axial mass flow is constant but the 
reactor walls are no longer cylindrical. The two additional limitations in applying the 
results of the calculations to gaseous nuclear reactor geometry are the use of the 
boundary layer equations close to the jet exit and the absence of an end wall in the cal- 
culations. Thus, the results of the calculations are approximations to reactor condi- 
tions. 

The amount of fuel present in a given volume can be compared with the minimum 
and maximum amounts that could be present. Fuel content would be a minimum if the 
jet and the coaxial stream were perfectly mixed before injection into the reactor. 

Thus, 


i rV 2 




If. 


r dr dx = 1 (py) 
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(39) 


If the jet and the coaxial stream are composed of pure component 1 and pure component 
2, respectively, 


so that 


Then, 
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and 


min 



(43) 


The amount of fuel would be a maximum if the initial jet were to remain a cylinder of 
uniform concentration with the same radius as the jet discharge tube. In this case 


max 




(44) 


Numerical Method 

In this section, a detailed analysis is presented of the numerical techniques used to 
solve the equations of coaxial turbulent jet mixing. Figure 2 is a general flow diagram 
for the numerical solution. The initial difficulty, from a numerical standpoint, is the 
application of the boundary conditions u = u and y = y as — °°. This difficulty was 
overcome by defining a parameter i such that as not only do the functions 

approach the boundary conditions, but also the derivatives of the functions are restricted 
to fall below some arbitrarily small parameter. Furthermore, since i//^ can be a 
numerically large value, a transformation was performed on the independent variable \]s 
to limit the range of integration from 0 to 1. An implicit finite -difference technique, 
the Crank-Nicholson method (ref. 9), was employed to solve the system of parabolic 
equations. Stability is inherent in such a scheme, and a high degree of accuracy can be 
obtained by a judicious choice of interval size. 

Let 




(45) 


so that 


j_ _ _i a 

^oo d: P 


(46) 


Applying this linear transformation to equations (20), (21), and the inverse of (25) results 
in the following equations: 
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( 47 ) 
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Correspondingly, the initial conditions become 


u = 1 0 < i// < — — 

2<K 


u = u ip > — 

2*. 


(48) 
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The boundary conditions are transformed to 


(51) 
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Finite-Difference Equations 


Consider a linear parabolic equation of the form 
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(54) 
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Equations (47) and (48) can be considered of this form if 

A (ffi.x) (55) 

U .,.2 


and 

Affix) = (56) 


and if the values of the quantities in equations (55) and (56) are initially assumed to be 
those at the previous axial position. By successively solving equation (54) and recomput- 
ing A u and Ay until no further change occurs, the correct values are obtained. 

Now consider a net R, ^ constructed on the region of interest. Let the subscript i 
represent discrete points on the (//-coordinate, and j discrete points on the x- 
coordinate. The points of the (//-coordinate are constructed such that i/^ = (i - l)Ai// 

with i ranging from 1 to N. Thus, the notation u. ■ corresponds to the functional 
r— . J 

value of u 0 /a, Xj). If forward differentiating is used over intervals in the x direction, 

then equation (54) can be integrated between mesh points to yield 



dl// = 




(57) 


If the integrand in equation (57) remains constant over the small interval Ai//, 


*£(u. . , - u. J=/a*A 
Ax 1,3+1 1,3 1 d\p).l • 1 

\ /1H — ,JH — 

2 2 


-A 


du\ 

d */U 

2 2 


(58) 


Central differences are used for the right side of equation (58); values of A and u at 
half intervals are approximated by the average over the whole interval; and the abbrevia- 
tions 
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( 59 ) 


w?= A * * 

1 4A \}y 


are introduced. Substituting these results in equation (58) yields 


'iVl, j+l -(*i * w i‘ + f^)“l,j + l + w t“i + l,j+l 


— - wTu. 1 . - l - wT - w1|u. • - wlu. 1 . (60) 
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Equation (60) describes a linear set of equations for the unknowns u. . Two 

J+ A 

additional equations are needed to complete the set in which i ranges from 1 to N. 

Integrating equation (54) from i = 1 to i = 1— and using the fact that (du/di//)-r-_ n = 0 

2 

gives 


( U 1 i+i " u i ))(—)= ( A ^=-) 

1,1+1 1,1 \2A xj y tyjl 1 
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Following the same procedure used to obtain equation (60), equation (61) can be ex- 
pressed as follows: 


-K + Kj+i + W 1 u 2, j+l = - U 2, j 


— - wj J “i j 

?At 1 


Applying the boundary condition at 4 / = 1 (i. e. , i = N) by using equation (60) with i = N 
and noting that u N+ j j + i = u e yields 

W N U N-1, j+l " ( W N + W N + “r) u N, j+l " 2w N u e " w n u N- 1, j "(t7 " W N _ w n) u N, j (63) 


The relations given by equations (60), (62), and (63) are written in matrix notation as 


follows: 


?Vi = d 


where the vectors u^ + ^ are the unknown quantities. The vector d has the following 
elements: 


II ii I I i 



The tridiagonal matrix B, whose superdiagonal and subdiagonal elements are wt and 
wr, respectively, has the following diagonal elements: 



The solution of equation (64) is directly obtained by Gaussian elimination of the subdiag- 
onal elements and a back substitution to obtain u j + j- 

Now that a method to solve equation (54) has been devised, the calculational pro- 
cedure for solving equations (48) and (49) is as follows: 

(1) Initially, at an axial length Xj, use values of r, p, u, and y from Xj -^ to 
compute A u and Ay. 

(2) Solve for Uj +1 and yj +1> 

(3) Recompute the coefficients A u and Ay. 

(4) Iterate between steps 2 and 3 until the change in A and A is less than a 

u y 

specified amount. 

This procedure was programmed for the IBM 7094 n 7044 Direct Couple System in 
FORTRAN IV. A listing of the program is given in appendix B. For most cases con- 

_ Q 

sidered, a value of At// = 1/200 and an initial Ax = 10“ gave satisfactory results. 
However, at larger axial positions, larger values of Ax can be used because of the 
decaying effects of the initial step profiles. An heuristic approach was used to alter Ax; 
if the iteration procedure converged in three or less iterations, Ax was increased by 
0. 5 percent but was limited to a value of 0. 4. Running time, of course, varied with 
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initial input; however, an average time of approximately 0. 3 minute per unit x could be 
expected. 


Program Input and Output 

The input required for a calculation consists of the following information: 

(1) The initial ratio of coaxial-stream velocity to jet velocity u 

V 

(2) The mass fraction of component 1 in the initial jet yj and in the initial coaxial 

stream y. 

J e 

(3) The ratio of the molecular weight of component 2 to the molecular weight of com- 

ponent 1 mg 

(4) The constants in the eddy viscosity formulation kj and kg 

(5) The turbulent Schmidt number Sc^ 

(6) The ratio of reactor diameter to jet radius d 

(7) The constants in the reference density formulation A and B 

(8) The axial positions at which output is desired x 

The output listing reproduces the input and, thereafter, at each axial position, gives 

(1) Axial position x 

(2) Eddy viscosity e 

(3) The product of density and eddy viscosity pe 

(4) The following information for axial velocity, mass fraction, and mole fraction: 

(Centerline value - Coaxial-stream value)/(l - Coaxial -stream value) 

The radial variations with stream function, also converted to radial position r, and the 
ratio of radial position to half radius r/r^yg for axial velocity, mass fraction, and 
mole fraction are provided in the following form: 

(Local value - Ambient stream value) 

(Centerline value - Ambient stream value) 

The following quantities are also listed: 

2 

(1) Momentum flux normalized with the square of the centerline velocity r/u£ 

(2) Mass flux normalized with the product of centerline velocity and mass fraction 

(3) Eddy viscosity and the product of density and eddy viscosity divided by their 
centerline values e/e^ and pe/(pe)^ 

Both sides of the centerline compatibility condition are printed next, followed by the 


16 


values of velocity and density at the largest stream function position in the calculation. 
Finally, the "line of sight" concentration c*, the dimensionless mass of component 1 I 
and the ratio of mass of component 1 to initial mass ft// ratio) are listed. 


RESULTS AND DISCUSSION 

Sample results from the computer program are presented to illustrate the mixing of 
a heavy, slow-moving jet and a lighter, faster-moving coaxial stream. First, however, 
a limiting case is discussed. 

Schlichting (ref. 7, p. 607) presents a similarity solution for an isothermal, tur- 
bulent free jet mixing with a quiescent ambient stream of the same composition that is 
valid far downstream. Schlichting' s value for the proportionality constant kg in the 
far -jet eddy viscosity formulation was adopted in order that the numerical solution repre 
sent this limiting case. The value of the constant of proportionality in the near -jet 
eddy viscosity formulation was chosen so that the numerical solution agreed with the 
similarity solution far downstream. 

Figure 3(a) shows a comparison of centerline velocity and half radius calculated 

from the similarity solution and the results of the numerical solution using two different 

-3 

values of kj. A value of k^ = 0. 75x10 leads to good agreement downstream and was 
therefore used in subsequent calculations. Radial velocity profiles rapidly change from 
the initial step profile and gradually merge into the similarity profile. Figure 3(b) 
illustrates that at an axial position of 16 jet radii the profile is almost similar, whereas 
at 50 jet radii the agreement is essentially exact. In these and subsequent calculations, 
the near jet formulation for eddy viscosity was used until (u^ - u e )/(l - u g ) = 0. 99; 
thereafter the far jet formulation was used. The calculation was repeated for 0.98 in 
order to see whether the particular choice of 0. 99 was critical. Although there was 
some difference initially, the difference between the two solutions was soon indistin- 
guishable. 

Choosing the proportionality constants in the eddy viscosity formulations by the 
method just discussed leads to a discontinuity in eddy viscosity at the axial position 
where one formulation replaces the other. It was felt that agreement with the similarity 
solution was more important than a continuous variation of eddy viscosity. Figures 4 
and 5 show typical calculations of the axial and radial variations of eddy viscosity and 
the product of density and eddy viscosity, respectively. In both cases, the radial 
variation is much greater in the near jet than farther downstream. The product of 
density and eddy viscosity varies less in the radial direction than does eddy viscosity. 

The effects of two different initial velocity ratios on velocity and mass fraction 
profiles are illustrated in figures 6 and 7. The higher velocity ratio, of course, re- 
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suits in a more rapid decay in centerline values and a narrower jet. In all calculations, 
the generally accepted value of 0. 7 (ref. 5, p. 422) was used for the turbulent Schmidt 
number. 

The effects of two different initial density ratios on velocity and mass fraction pro- 
files are illustrated in figures 8 and 9. As expected, the lighter jet decays faster and 
results in a narrower jet. 

The influence of reference density on velocity and mass fraction profiles are shown 
in figures 10 and 11. The formulation based on centerline density leads to much more 
rapid decay of centerline values and to a narrower jet than does the formulation based on 
ambient stream density. The formulation based equally on centerline and ambient 
stream densities falls between these results. These figures demonstrate that it will be 
necessary to determine experimentally if any of these formulations are adequate. 

Figures 12 and 13 show the effect of initial velocity ratio and initial density ratio 
on the mass of the major jet component contained within a given volume. The maximum 
and minimum values are included for comparison. 


CONCLUDING REMARKS 

A computer program to describe isothermal, turbulent jet mixing of two gases was 
written using the axisymmetric forms of the boundary layer momentum and diffusion 
equations. The coaxial stream is considered to be infinite in extent. Eddy viscosity is 
represented by an expression that provides for both radial and axial variation. Typical 
running time is less than 1 minute to produce time -averaged velocity and mass fraction 
distributions. 

Experimental data are required for further progress. Constant-density experiments 
at large initial velocity ratios will determine if the numerical value of kg used is 
appropriate. If not, the large initial velocity ratio data can be used to determine a new 
value. The value of kj can also be obtained from the same experiments. With the 
values of kj and kg determined, variable -density experiments at large velocity ratios 
can be used to determine a suitable reference density in the eddy viscosity formulations. 

In the comparison of computer results and experimental jet-mixing data, the initial 
conditions on the equations may have to be modified to account for the finite wall thick- 
ness of the jet discharge tube and for the distribution of velocities in the initial jet and- 
the coaxial stream. For gas -core nuclear reactor calculations, the absence of an end- 
wall boundary condition is probably a serious restriction. In addition, the assumption of 
constant pressure and the use of the boundary layer equations near the jet exit are 
approximations. These restrictions can be removed by using the full Navier -Stokes 
equations rather than the boundary layer equations. However, the problem of turbulence 
and a method to characterize an eddy viscosity remain. 
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For unbounded turbulent jet mixing, the computer program provides a rapid solution 
that reduces to the similarity solution when the density is constant and the coaxial 
stream is quiescent. Improved values for the proportionality constants in the eddy 
viscosity formulations, or an entirely new expression, can easily be incorporated. 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, September 8, 1967, 
122-28-02-16-22. 
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APPENDIX A 


V A y 


SYMBOLS 

constants in reference density formula 
quantities defined by eqs. (55) and (56) 
coefficient matrix 
elements of B matrix 
mole fraction of component 1 

dimensionless mole fraction of component 1, C/Cj 
line of sight concentration 
reactor diameter 

dimensionless reactor diameter, D/R^; elements of d vector 
vector defined by eqs. (65) 

rl rd/2 

dimensionless mass of component 1,1 / pyr dr dx 

Jo Jo 

constants in eddy viscosity formulations 
reactor length 

dimensionless reactor length, L/Rj 
molecular weight 

dimensionless molecular weight, M/Mj 
upper limit of i 

radial position; finite difference net 
dimensionless radial position, R/R.. 

dimensionless half radius (i. e. , position at which (u - u )/(l - u ) = 1/2) 
turbulent Schmidt number 
axial velocity 

dimensionless axial velocity, U/Uj 
velocity vector 
radial velocity 
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v dimensionless radial velocity, V/Uj 

3 

W amount of fuel contained within a given volume, 27rPjYjRj I 

w abbreviation defined by eq. (59) 

X axial position 

x dimensionless axial position, X/Rj 

Y mass fraction of component 1 

y dimensionless mass fraction of component 1, Y/Yj 

E eddy viscosity 

e dimensionless eddy viscosity, E/UR- 

1 J 

e* dimensionless constant-density eddy viscosity 

p density 

p dimensionless density, P/Pj 

x [/ stream function 

Ip normalized stream function, xfz/xp^ 

xf/^ maximum value of xjy 

Subscripts: 
av average 

<t centerline 

e coaxial stream 

i point on xf/ coordinate 

j jet; point on x coordinate for numerical solution 

max maximum 

min minimum 

0 reference 

1 major component of jet 

2 major component of coaxial stream 
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YED1581 C A TODD COAXIAL FLOW 

A COMPUTER PROGRAM FOR CALCULATING ISOTHERMAL 
TURBULENT JET MIXING OF TWO GASES 

SYMBOL DEFINITION 

AP THE A COEFF. FOR THE DIFFUSION EON. 

AU THE A COEFF. FOR THE MOMENTUM EON. 

BP THE B COEFF. FOR THE DIFFUSION EON 

BU THE B COEFF. FOR THE MOMENTUM EON. 

DPSI THE INTERVAL SIZE IN THE PSI DIRECTION 

DX THE INTERVAL SIZE IN THE X DIRECTION 

EA COEFF. A , FOR REFERENCE DENSITY CALCULATION 

EB COEFF. B r FOR REFERENCE DENSITY CALCULATION 

EPS RHO TIMES THE EDDY VISCOSITY 

I INDEX VARIABLE 

I HALF VALUE OF I WHEN R-R-l/2 
ITER THE ITERATION COUNTER 
J INDEX VARIABLE 

K A COUNTER TO CONTROL THE OUTPUT 

KGO A LOGICAL VARIABLE TO CONTROL THE CALCULATION OF THE 
CONTAINMENT FACTOR 

NCOPY THE NUMBERS OF COPIES WANTED — NORMALLY l 

NPTS THE NUMBER OF POINTS IN THE PSI DIRECTION 

NPTSX THE NUMBER OF X VALUES TO BE OUTPUTED 

NORMPO THE EUCLIDIAN NORM OF RHO AT X-DX 

NORMP2 THE EUCLIDIAN NORM OF RHO AT X 

NORMUO THE EUCLIDIAN NORM OF U AT X-DX 

NORMU2 THE EUCLIDIAN NORM OF U AT X 

PMAX PSI-INFINITY AT X 

PMAXOD PSI-INFINITY AT X-DX 

PSI THE INDEPENDENT VARIABLE 

PSIMO PSI-1/2 AT X-DX 

PS I Ml PSI-1/2 AT X 

PTS THE NUMBER OF INTERVALS IN THE PSI DIRECTION 

PTSX THE NUMBER OF OUTPUT X'S 

PZERO THE PSI-RATIO 

RHALFO R-l/2 AT X-DX 

RHOE THE EDGE DENSITY 

RHOR REFERENCE DENSITY 

RHOO VALUE OF RHO AT X-DX 

RHOl A GUESSED VALUE OF RHO AT X 

RH02 A COMPUTED VALUE OF RHO AT X 

RO VALUE OF R AT X-DX 

R 1 A GUESSED VALUE OF R AT X 

R2 A COMPUTED VALUE OF R AT X 

SC THE TURBULENT SCHMIDT NUMBER 

SD A DIAMETER 

SM2 A MOLECULAR WEIGHT RATIO 

SUM TEMPORARY STORAGE 

SUM I TEMPORARY STORAGE 

TOL THE TOLERANCE TO TERMINATE THE ITERATION 

UE THE EDGE VELOCITY 

UO U AT X-DX 

UM INTERPOLATED U AT X-DX 

U1 A GUESSED VALUE OF U AT X 

U2 A COMPUTED VALUE OF U AT X 

XI THE CONTAINMENT FACTOR 

XO THE PREVIOUS VALUE OF THE INDEPENDENT VARIABLE, X 

XI THE CURRENT VALUE OF X 
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XK1 A TURBULENCE FACTOR FOR EDDY VISCOSITY IN THE NEAR JET 

XK2 A TURBULENCE FACTOR FOR EDDY VISCOSITY IN THE FAR JET 

XM A CONTROL VARIABLE 

XM=0 • NORMAL INPUT 

XM= 1 • NORMAL INPUT PLUS RESTART DUMP CARDS 
XPCH VALUES OF X TO BE OUTPUTED 
YE THE EDGE VALUE OF Y 

YJ THE JET VALUE OF Y 

YM INTERPOLATED Y AT X-DX 

YO VALUE OF Y AT X-DX 

Y 1 GUESSED Y AT X 

Y2 COMPUTED Y AT X 

INPUT 


UE 

CARD 

YJ 

1 FORMAT 

YE 

8F10.X 

SM2 

XK 1 

XK2 

SC 

TOL 

PTS 

CARD 

PTS X 

2 FORMAT 

PMAX 

8F10.X 

SD 

EA 

EB 

XM 

NCOPY 


CARD 

3+ FORMAT 

8FI0.X 





XPCH(l) 

XPCH ( 2 ) 

— 

ETC 

— 



XPCH ( 8 ) 


COMMON/ DATUM/ RO*UQ*RHOG*PSI* NPTS , UE , Y J , YE r SM2 , XK I , 

1XK2,SC,RH0E ? DPSI , YO, RH02 , U2 ♦ R2 , PS I Ml * RHOE 1 ,RHALF1,X1, IHALF,Y2 
COMMON RHOR,EPS( 500) 

DIMENSION U0( 500) , U1 ( 500 ) , U2 ( 500 ) , RHOO ( 500 ) , RHO 1 ( 500 ) , RH02 ( 500 ) , 
IRC) (500) r Rl (500) ,R2 ( 500) ,VL( 500) , V ( 500 ) , BU ( 500 ) , 

2BP( 500) ,XPCH( 50) ,PSIPCH< 50) , PS I ( 500) 

DIMENSION YO ( 500 ) ,Y1(500) ,Y2(500) 

DIMENSION YM(500) ,UM(500) 

REAL N0RMU0,N0RMU1 ,NORMU2,NORMPO,NORMP1 ,N0RMP2 
SUB(XltYl,X2,Y2,X)=(X-Xl)/(X2-Xl)*(Y2-Yl)+Yl 
COMMON XO ( I ) t PMAX ,NORMUO , NORMPO , DX , XZ t RHALFO , X I , K ,PZERO 

READ ALL REQUIRED INPUT DATA AND INITIALIZE VALUES OF 

STARTING DX. 


1 READ(5,400)UE t YJ,YEtSM2tXKl,XK2tSC f T0L,PTS,PTSX ,PMAX ,SD ,EA,EB 
1 9 XM t TIB 
NCOP Y = T I B 
M = XM 

PMAXOD=PMAX 

NPTS=PTS+I. 

NPTSX=PTSX 
T=SM2— 1 • 

RH0E=(YJ*T+1. )/ ( YE*T+i.) 

DPS I = 1 • / PTS 
K = I 

READ ( 5 * 400 ) ( XPCH ( I ) t I =1 , NPTSX ) 

X I =0 * 

DO 32 I = 1 » NCOP Y 
32 CALL INITAL(M) 

KGO=I 

IF(SM2.GT.l.) KG0=2 
WRITE(6,471) SD ,EA t EB 
IF(M.NE.O) GO TO 100 
DO 9969 I = 1 ♦ NPTS 
J = I 

9969 I F ( RO ( I ) .GE.SD/2. ) GO TO 9970 

9970 PZ ERO = PS I ( J ) 
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I F ( J • EQ • NPTS ) PZER0=PZER0+RHUE*UE/2.*< SD**2/4.-R0 ( NPTS ) **2 > 
1/PMAX 

471 F0RMAT(2HKD,G13.5,2X,1HA,G13«5,2X,1HB,G13.5 ) 

400 FORMAT { 8F10*7 ) 

100 NORMP 1=N0RMP0 
N0RMU1 =N0RMUQ 
PMAX=PMAX0D 

SET FIRST GUESS OF ( X+DX ) =U < X ), Y ( X+DX ) =Y < X ), AND R(X+DX)=R(X) 

DO 10 I = 1 1 NPTS 
RH01U ) = RHOO ( I > 

U1 ( I ) =UO ( I ) 

Y 1 ( I ) = YO ( I ) 

10 R 1 ( I ) =RO ( I ) 

X 1 =X0+DX 
I TER = 0 
NPTS 1 =NPTS 

IF PSI-MAX HAS CHANGED INTERPOLATE VECTORS TO CORRESPOND 

TO NEW LENGTH* 

101 DO 102 1=1, NPTS 
YM< I ) = YO ( I ) 

102 UM< I )=U0( I ) 

IF(PMAX.EQ.PMAXOD) GO TO 50 
DO 103 1=1, NPTS 

VI ( I ) = PS I ( I >*PMAX 

103 V( I ) =PS I ( I )*PMAXOD 
DO 104 I = 1 , NPT S 

CALL SINTP(V,U0,NPTS,V1 ( I ) ,UM( I) ) 

104 CALL SINTPt V,Y0,NPTS,V1 ( I ) ,YM( I ) ) 

COMPUTE PSI-MAX AND PSI-1/2 AND RHO-EPSILON 

50 TEST=(U1 ( 1 )— UE ) / C 1 • — UE ) 

RH0R=EA*RH01(1)+EB*RH0E 
IF(TEST.LE*.99)G0 TO 12 
RH0E1=XK1*X1*ABS<UE-1. ) 

I HAL F = 0 

DO 90 1 = 1, NPTS 

90 E PS ( I ) = ( RHOR/ RH01 ( I ) )**2*RH0E1 
GO TO 61 

12 TEST=.5*(UE+U1(1) ) 

IFtUE.GT.l.) GO TO 490 
DO 14 1=1, NPTS 
I HAL F = I 

I F { U1 ( I UGE.TEST) GO TO 14 
GO TO 15 
14 CONTINUE 

490 DO 491 1=1, NPTS 

I HAL F= I 

I F ( U 1 ( I ) *LE . TEST ) GO TO 491 
GO TO 15 

491 CONTINUE 


COMPUTE VALUES OF AU,AP,BU,AND BP AND SOLVE FOR NEXT 
APPROXIMATION OF U,Y ,RHO,AND R. 


15 DO 456 1 = 1,1 HAL F 
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456 V( I ) =PMAX/U1 { I ) 

CALL FNTGRH IHALF*DPSI,V»V1 ) 

RHALFl=SQRT(2./RHOR *Vi ( I HALF) ) 

RH0E1=XK2*RHALF1 *ABS(UE-U1 ( 1 ) ) 

13 DO 17 I = 1 » NPTS 

17 V ( I ) = PMAX / U 1 ( I) 

CALL FNTGRLtNPTSfDPSI ,V*V1) 

EPS( 1 >=RH0E1*RH0R/RH01 { 1 ) 

DO 62 1=2, NPTS 

62 EPS ( I )=2.*RH0R*RH0E1/RH01 ( I )**2/Rl ( I )**2*V1 ( I ) 

61 DO 455 1=1, NPTS 
AU = 1. 

A P= 1 • 

BU( I ) =EPS ( 1 )*RH01 ( I ) *#2*U1 ( I ) *R 1 ( I >*#2 
BU ( I ) =BU ( I ) /PMAX**2 
455 B P ( I )=BU( I )/SC 

CALL SOLVE ( AU , BU , NPTS , UE , U2 , NORMU2 , DX , DP S I ,UM) 

CALL SOLVE ( AP ,BP , NPTS ,YE , Y2 ,NORMP2 ,DX,DPSI , YM > 

T=SM2— 1 • 

T1=1./YJ 

DO 7777 1=1, NPTS 

7777 RH02( I )=(T+T1)/(Y2( I ) #T +T1) 

NPTS2 = NPT S 1 
DO 18 1=1 ,NPTS2 

18 V( I )=PMAX/RH02( I )/U2( I ) 

CALL FNTGRL (NPTS2,DPSI , V , VI ) 

DO 19 1=1, NPTS2 

19 R2U )=SQRT(2. *ViU ) ) 

DO 11 1=1, NPTS 1 

11 V ( I ) = ( U1 { I ) -UE ) / ( 1 .-UE) 

CALL FMTGRL(NPTS1,DPSI,V,V1) 

DO 7070 1=1, NPTS 
11 = 1 

I F { V 1 ( I >.GT..495*PMAX) GO TO 7071 

7070 CONTINUE 

7071 PS I Ml = PS I (ID 
SUM I =0 * 

DO 460 I =2 , NPT S 
GO TO (461,462) ,KGO 

461 SUM=Y2 ( I — 1 ) 

GO TO 463 

462 SUM= 1 • — Y 2 ( 1-1 ) 

463 IF(R2< I ) .GT.SD/2. ) GO TO 460 
SUMI=$UMI+SUM*RH02 ( I -1 ) *R2 ( 1-1 ) * ( R2 ( I )-R2 ( 1-1 ) ) 

460 CONTINUE 

CHECK TO SEE IF CONVERGENCE CRITERIA HAS BEEN MET. 

I F ( ABS { ( N0RMU2— N0RMU1 ) /N0RMU2 ) • GT • TOL ) GO TO 20 
I F ( ABS ( ( N0RMP2— NORMP 1 ) /N0RMP2 ).GT.TOL) GO TO 20 
TEST=(U2 (NPTS) -U2( NPTS- 1 ) )/DPSI 
IF(TEST.GT,.001*PMAX) GO TO 70 
DEBUG XI, DX, PMAX 
DEBUG IPSI ( I ) ,I=lrNPTS,20) 

DEBUG <R2 (I), 1=1, NPTS, 20) 

DEBUG ( U2 ( I ) , 1=1 , NPTS, 20) 

CONVERGENCE CRITERIA HAS BEEN MET 
CHECK WHETHER DX CAN BE INCREASED 
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CALL TIMLFT(TIMI) 

C 

C IF THE TIME REMAINING IS LESS THAN .1 MIN* OUMP FOR RESTART. 
C 

IF!TIMl/3600..GT..l) GO TO 3333 
XZ=K 

CALL BCDUMP ( XQ! 1 ) * XO ( 8 ) ) 

CALL BCDUMPtPSI ( 1) * PS I (NPTS) ) 

CALL 8CDUMP(R0f 1 ) f RO(NPTS) ) 

CALL BCDUMP ( UO ( 1 ) * UO ( NPTS ) ) 

CALL BCDUMP ( RHOO ( 1 ) * RHOO (NPTS) ) 

CALL BCDUMP! YOU) t YO(NPTS) ) 

STOP 

C 

C IF ITER IS LESS THAN 3, INCREASE DX. 

C 

3333 IF! ITER.LT. 3) DX=1.05*DX 
I F( DX.GT • .4) DX=.4 
X I = X I +DX*SUM I 

C SHOULD WE PUNCH OUT AT THIS X 

C 

IF(X1.GT.XPCH(K) ) GO TO 30 
C 

60 X0=X1 

PMAXOD=PMAX 
N OR MU 0= NOR MU 2 
NORMPO=NORMP2 
RHALFO=RHALFI 
PS I MO=PS I Ml 
RH0E0=RH0E1 
DO 21 I = 1 * NPTS 
UD( I ) =U2 ( I ) 

RHOO! I ) =RH02 ( I ) 

YO ( I ) =Y2 ! I ) 

21 RO ( I > =R2 l I ) 

GO TO 100 

C 

C NO CONVERGENCE 

C 

20 M0RMU1=N0RMU2 
N0RMP1=N0RMP2 
ITER=ITER+1 
DO 22 1 = 1* NPT S 
U1 ( I )=U2( I ) 

RH01 ( I )=RH02 ( I I 
Y 1 ( I ) =Y2 ( I ) 

22 R 1 C I ) =R2 ( I J 
GO TO 50 

C 

C PUNCH OUTPUT 

C 

30 DO 31 I = 1 » NCOPY 

31 CALL OUTPUT ( XPCH ( K ) ) 

WR I T E ( 6 * A 7 0 ) X I 

DO 9971 1=1, NPTS 

J = I 

9971 IF(R2( D.GE.SD/2.) GO TO 9972 

9972 TEST=PSI!J> 

IF! J.EQ.NPTS)TEST=TEST+RH0E*UE/2.*(SD**2/4.-R2 (NPTS) ^*2 ) 
1/PMAX 
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RATP=TEST/PZERO 
WRITE(.6,479) RATP 
WRITE (6,482) SUM I 
482 FORMAT { 5H I ( X ) , 2 X , G 1 3 . 5 ) 

479 FORMAT { 1 OHK PS I -RAT I 0 , 2 X , G 1 1 .4 ) 

470 FORMAT (2HKI , G1 3 # 5 ) 

K = K+ 1 

1FIX1 .LT.XPCH(K) )G0 TO 60 

GO TO 1 

MUST CHANGE PSI-INF 

70 PMA X = PMAX+2 » 

GO TO 101 
END 

BFTC SOLVE 

A ROUTINE TO SOLVE A PARABOLIC EQUATION BY THE CR AtyK-N I CHQL SONS 

ALGORITHM. 

SUBROUTINE SOL V E ( A , B , N , HMAX , H , NORMH , DX ,DPS I ,H0) 

REAL NORMH 

DIMENSION B ( 500 ) ,H(500) , SB ( 500 ) , SD ( 500 > , HO ( 500 ) ,SA(500> 

ISC (500) , WP ( 500 ) , WM { 500 ) 

N 1 =N— 1 
S A ( 1 ) =0 . 

T =A/4 . / DPS I 
T 1 = DPS I / DX 
DO 5 1=2, N 

5 WM(I)=T*(B(I)+B(I-1) ) 

DO 6 1=1, N1 

6 WP< I )=T*(B( 1 + 1 )+B( I ) ) 

SB ( 1 ) =— ( WP ( 1 )+T 1/2 • ) 

SC( 1 )=WP( 1 ) 

SD ( 1 ) = { WP ( 1 } —T 1 / 2 . ) *H0 ( 1 ) — W P ( i)*H0(2) 

N=N— 1 

S A ( N ) =WM ( N ) 

SB(N)=— (WP(N)+WM(N)+T1 ) 

SDlN)=-l.*WP{N)*HMAX-WM(N)*HO(N-l ) -( -WM ( N > + T 1 -W P ( N> )*HO(N) 
1-WP(N)*H0(N+1 ) 

N 1 =N— 1 
DO 7 1=2, N1 
SA( I )=WM( I ) 

SC( I ) =WP ( I ) 

SB ( I )=-( WM( I )+WP( I ) +T 1 ) 

7 SD ( I ) =— WM ( I )*H0( I-1)-WP( I )*H0( 1 + 1 >-(Tl-WM( I )-WP( I) )*H0( I ) 

DO 2 1=2, N 

SB( I )=SB( I )-SA( I )/SB( 1-1 ) *SC( I-l ) 

2 SD ( I ) =SD ( I ) -S A ( I )*SD( 1-1 )/SB( 1-1 ) 

H(N)=SD(N)/SB(N) 

N0RMH=H(N)**2 
DO 3 1=1, N1 
J=N — I 

H( J ) = ( SD( J )-SC ( J )*H( J + l ) ) /SB ( J ) 

3 NORMH=NORMH+H< J )**2 
N = N+ 1 
H ( N ) =HMA X 
NORMH=SGRT(NORMH) 

RETURN 
END 

S I BFTC OUTPUT 


SUBROUTINE OUTPUT(XX) 

COMMON/ DATUM/ RO ,UO , PO, PS I ,N1 , UE ,Y J ,YE ♦ SM2 r XK1 , XK2 , SC , RHOE , DPS I , 
lYO t P2,U2,R2,PMl,PEl T RH2,Xl, IH,Y2 
COMMON RHOR,EPS{ 500) 

DIMENSION P2 ( 500 ) , U2 ( 500 ) , R2 ( 500 ) ,PO(500) , UO ( 500 ) , RO ( 500 ) , 

IPS I ( 500) ,PSIP ( 500) t P X ( 500 ) ,UX ( 500) ,RX( 500) ,ROX< 500) ,URAT{500) , 
2YRATI 500) , POP ( 500) ,TAU(500) , XMU( 500) , Y0( 500) , YX( 500) , Y2( 500) 
DIMENSION ROXI ( 500 ) t EP( 500 ) ,RHOEP( 500 ) 

COMMON Zl t PMO 

SUB (XI ,Y1 ,X2,Y2,X)=(X-X1 ) / (X2-XL )*( Y2-Y1 )+Yl 
PMA=SUB(X0,PM0,X1,PM1,XX) 

DO 91 1=1 t N1 

91 PS I P ( I ) = PSI ( I ) * PMO 
DPSI=PSIP(2)-PSJP(1> 

DO 51 1=1, Nl 
IH = I 

T=(U2( I ) -UE ) / (U2( 1 ) — UE ) 

I F ( T *GT * *5 ) GO TO 51 
GO TO 52 

51 CONTINUE 

52 RH 1 =R2 ( I H ) 

50 PMX=1. 

S0T=1 • 

11 RHOEX=PE 1 /SOT 

RHOX = SUB( XO, P0( 1 ) , XI » P2 ( 1) ,XX) 

I F ( RHOE *NE • 1 • ) 

1RH0X=( RHOX-RHOE )/ (1 .-RHOE) 

UOX=SUB( XO,UO( 1 ),X1,U2(1),XX) 

URATX= ( UOX-UE ) / ( 1 .-UE) 

T =0 • 

IF ( SM2.NE. i • ) T=1./YJ/(SM2-1. ) 

YOX=SUB ( XO,YO ( 1 ) ,X1,Y2( 1 ) ,XX ) 

YR ATX=0. 

IFlYE.NE.l.) YRATX=( YOX-YE>/( l.-YE) 

W = S M 2 - 1 . 

W1=1./YJ 
DO 1 1=1 , N 1 

YXt I ) = SUB ( XD , YO ( I ) ,Xl r Y2( I ) ,XX) 

PX { I ) = (W+W1 ) / ( YXU )*W + W1 ) 

UX ( I ) =SUB ( XO,UO( I ) ,X1,U2 ( I ) ,XX ) 

R X ( I )=SUB( XO y RO( I ) ,X1 ,R2 ( I ) ,XX) 

ROX ( I ) =0 . 

I F ( IH.NE.O) 

1 R OX ( I ) = R X ( I ) / RH 1 
URA T ( I ) = ( UX ( I )-UE) / (UOX-UE ) 

XMU ( I ) = YX ( I )*PX ( I )*UX( I )**2 
POP( I ) =0 . 

IF( RHOE.NE.l. ) 

1 P 0P< I ) = ( PX < I J-RHOE )/ ( PX ( 1 ) -RHOE ) 

YR AT ( I ) =0 . 

IF ( YE.NE.YJ J YRAT ( I ) = ( YX( I ) -YE ) / ( YOX-YE ) 

ROXI ( I ) =RX( I)/XX 
1 CONTINUE 

DO 521 1=1, N1 

EP ( I ) =E PS ( I ) /EPS ( 1 ) 

521 RHOE P ( I )=PX< I )*EPS( I )/PX( 1 )/EPS( 1 ) 

TAU( 1 )=PMX*(UX(2)-UX< 1 ) )/DPSI 
T AU ( N 1 )=PMX*(UX ( N 1 )-UX (Nl-1 ) ) /DPS I 
XMU ( 1 ) = ( YX ( 2 )-YX( 1 ) ) /DPS I 



XMU (N1 ) = ( YX ( N 1 )— YX(N1—1 ) )/DPSI 

N3=N1-1 

00 2 1=2, N3 

TAUU )=PMX*(UX( I + 1)-UX( 1-1 ) ) /2./DPSI 

2 XMU ( I ) = (YX( 1 + 1 ) — Y X < 1-1 ) )/2./DPSI 
00 3 I =1 ,N1 

T =— RHGEX#PX ( I )*UX( I )*RX( I ) /SOT 

TAU ( I )=TAU( I )*T /U0X**2 

3 XMU( I )=XMU( I )*T/SC /YOX/UOX 
RHOEPX=PX( 1 ) *EPS ( 1 ) 

RHX=RH1 

WRITE (6,500) XX,URATX,YRATX,RHOX,EPS( 1 ) ,RHOEPX ,RHX 
I MOD= l 0 
K = 1 

TEST=PS I HF 
I STRT= 1 

32 DO 6 I=ISTRT,N1 
N = I 

I F < PS I P < I ) .GT. TEST ) GO TO 7 

6 CONTINUE 

7 MOD=N/ I MOD 
IF(MOD.LT.l) MOD= 1 
I F ( K . E0« 2 ) N=N+2 
DO 4 I=ISTRT ,N, MOD 
R XXX = RX ( I ) / SOT 
IF{SM2.EQ.1.)XMU< I )=0. 

T5 = Y X ( I )*PX< I ) 

4 WRITE(6,501) PSIP(I) ,RXXX ,ROX ( I ) ,URAT( I ) , YRAT ( I ) ,POP( I ) , T AU ( I ) 

1 , XMU ( I ) , EP ( I ) , RHOEP ( I ) ,T5 

GO TO (30,31) ,K 

30 I STRT = N+ 1 
TEST=PMA 

I MOD= 1 0 
K = 2 

GO TO 32 

31 D1 = U0X*(U2( 1>-U0( 1 ) >/ (Xl-XO) *PX(1> 

D2=A./RX(2)**2*RH0EX*(UX(2)-UX( I) ) *RHOR 
T = ( UOX— UE ) / ( 1 • — UE ) 

IF(T.GT..99) D 1 =0 1 *PX ( 1 ) 

I F ( T .GT . .99 ) D2 =D2*RH0R 
WR I TE ( 6 , 600 ) D 1 , D2 

600 FORMAT ( 32HKCENTERL INE COMPATIBILITY VALUES 2G15.5 ) 

WRITE! 6,601 ) UX ( N 1 ) , P X ( N 1 ) 

601 FORMAT ( 7H UMAX= G 1 3 . 5 , 3 X , 9HRH0-M A X = G13.5 ) 

N 3 = N 1 — 1 
C STAR = 0 • 

I F ( RHOE • EO. 1 • ) GO TO 602 
DO 66 1=1, N3 
DR = R X ( 1 + 1 ) — RX ( I > 

66 CSTAR=CSTAR+(PX( I ) -RHOE )/( 1 .-RHOE ) *DR 

602 CONTINUE 

WRIT E(6, 520) C STAR 
520 F0RMAT(3HKC*,G13.5) 

K = 1 

RETURN 

500 FORMAT (19H1 AX I AL-L ENGTH , X , G1 1 .4 , 5X , 14H ( UO-UE ) / ( 1-UE ) ,G11.4, 

15X, 14H(Y0-YE)/( 1-YE ) , G 1 1 . 4 , 14X , 14H < PO-PE ) / ( 1-PE),G11.4 / 
213X,6HEPS-0 ,G11.4,12X,7HRH0EPS0,G11.4,14X,5HR-1/2,G11.4,12X, 
3// 

4 23H U-RATI0=(U-UE)/ (UO-UE)/ 23H Y-R AT I 0= ( Y-Y E ) / ( YO-Y E ) / 
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S23H P — RAT 10= ( P — PE ) / ( P0— P E ) / 

6 , 26H TAU NORMALIZED BY 1/U0**2 / 25H MU NORMALIZED BY 1/UO/YO // 
54H PS I ,7X,1HR, 10X,9HR/ (R-l/2) ,2X, 7HU-RAT I 0, 4X, 7HY-RAT I0,4X , 
69HP-RATI0 ,2X,3HTAU,8X,2HMU,9X, 

7 12HEPS/ ( EPS— 0 ) T lX f 1 7HRH0E PS/ ( RHOEPS— 0 ) , 

22 X r 5HRH0^Y ) 

501 FORMAT ( 10G11.4,6X,G11*4) 

END 

$ I B FTC INITAL 

SUBROUTINE INITAL (M) 

COMMON/ DATUM/ R,U,RHO,PSI , N » UE , Y J , Y E , SM2 t XK 1 , XK2 , SC , RHOE , 

1 DPS I f YO 

DIMENSION R( 500) ,U(500> ,RHO< 500) ,PSI (500) , V ( 500 ) , V 1 < 500 ) ,Y0( 500) 

REAL NU ,NP 

COMMON RHOR , EPS ( 500 ) 

COMMON X(n*PMO,NU,NP,DX,XK,RHFO,XI t K f PZERO 
IFIM.EO.l ) GO TO 100 
C 

T= ( SM2— 1 • ) #YJ 
PS I ( 1 >=0. 

NU = 0. 

NP =0 • 

DO 1 I =2 , N 

1 PSI( I )=PSI ( 1-1 )+DPSI 
DO 8 1 = 1, N 

IF(PSI(I).LE..5/PM0) GO TO 4 
U( I ) =UE 
RHO ( I >=RHOE 
GO TO 3 

4 U ( I ) =1 • 

RHO ( I ) = 1 * 

3 NU=NU+U( I )**2 

yo( r ) = i • 

IFIT.NE.O. ) YO ( I ) = { Y J/RHO ( I )#(!•+ l./T) — YJ/ T ) / Y J 
NP=NP+YO ( I )**2 
8 CONTINUE 
NU=SGRT(NU) 

NP=SORT ( NP ) 

20 CONTINUE 
51 DO 5 1=1, N 

5 V ( I )=2./RH0( I )/U( I ) *PMO 
CALL FNTGRL(N,DPSI,V,V1) 

DO 7 I = 1 , N 
7 R( I )=SORT( VI ( I ) ) 

RHF0=0. 

I F ( M . EG • 2 ) RETURN 
X = 0 

DX=1 • E— 2 
GO TO 101 

100 CALL BCREA0(X(1) ,X(9) ) 

CALL BCREAD(PSI (1 ),PSI (N) ) 

CALL BCREAD(R( 1 ) ,R(N) ) 

CALL BCREAD(U(1 ) ,U(N) ) 

CALL BCREAD(RHO( 1 ) ,RHO(N) ) 

CALL BCREAD(Y0(1 ) ,YO(N) ) 

PMO=X ( 2 ) 

NU=X< 3) 

NP=X ( 4 ) 

DX=X( 5) 

K =X ( 6 ) 



I 

1 


RHFO=X ( 7 ) 

X I = X ( 8 ) 

PZ E RO=X ( 9 ) 

DEBUG (X(I),I=1,8) 

DEBUG ( YO ( I ) , I =1 , N ) 

101 WRITE (6,500)UE,YJ,YE,SM2,XK1,XK2,SC ,X 

$OT = SQR T ( PMQ ) 

M0D=N/20 

DO 6 1=1, N, MOD 

PS I P = PS 1(1) *PMO 

RP = R( I ) 

6 WRITE(6,501 )PSIP,RP,U( I ) ,RHO( I ) 

RETURN 

500 FORMAT ( 1H1 ,30X,30HINPUT FOR TURBULENT JET MIXING // 
17X,2HUE,G11.4,7X,2HYJ,G11.4,7X,2HYE,G11.4,7X,2HM2,G11.4,7X,2HK1, 
2G11 .4,7X,2HK2,G11 .4/ 7 X , 2 HSC ,G 1 1 .4 // 

330X,24HINITIAL PROFILES FOR X= G11.4 // 

54H PSI ,9X,1HR,12X,1HU,12X,3HRH0 ) 

501 FORMAT ( 4G13.5) 

END 

$ I BFTC SINTP 

SUBROUTINE S INTP ( X ,Y ,N , XI , Y1 ) 

DIMENSION X ( 500 ) , Y ( 500 ) 

C 

DO 1 1=1, N 

K = I 

IF ( X1.GT.X( I )) GO TO 1 
IF (Xl.EO.Xl I) ) GO TO 2 
IF (Xl.LT.X(D) GO TO 3 

1 CONTINUE 

2 Y 1 = Y { K ) 

GO TO 100 

3 IF (K.EQ.l) K=2 
IF IK.EO.N) K=N— 1 

I F ( Y ( K — 1 ) .NE.Y(K) ) GO TO 5 
Y 1 = Y ( K ) 

RETURN 
5 CONTINUE 

W1 = (Xl-X(K)) *(X1-X(K+1) ) / « X C K— 1 )—X ( K ) ) / ( X < K- 1 )-X ( K+ 1 ) ) 

W2 = <Xl-X(K-l))*<Xl-X(K+l>)/(X(K)-X(K-l)>/(X(K)-X(K+m 
W3 = ( X 1— X ( K— 1 ) ) * ( X 1—X ( K ) )/ ( X(K+1 ) — X (K — 1 1 )/ ( X ( K+ 1 ) — X { K ) ) 

Y 1 = Y ( K — 1 )^W1 + Y(K)^W2 + Y(K + 1 )*W3 
100 RETURN 
END 
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Axial position, x 

(a) Centerline velocity and half radius. 



(b) Velocity profile. 


Figure 3. - Comparison of numerical and similarity 
solutions. Quiescent coaxial stream; constant 
density. 
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0 2 4 6 0 .4 .8 1.2 1.6 


Axial position, x Radial position, r 

(a) Axial variation. (b) Radial variation. 

Figure 4. - Typical variation of eddy viscosity. Velocity ratio, 30, density 
ratio, 0.170, reference density, average of centerline density and 
coaxial-stream density. 



0 2 4 6 0 .4 .8 1.2 

Axial position, x Radial position, r 


(a) Axial variation. (b) Radial variation. 

Figure 5. - Typical variation of product of density and eddy vis- 
cosity. Velxity ratio, 30, density ratio, 0. 170, reference 
density, average of centerline density and coaxial -stream 
density. 
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0 .4 .8 1.2 0 .4 .8 1.2 

Radial position, r 

(a) Velocity ratio, 10. (b) Velocity ratio, 30. 

Figure 6. - Effect of initial velxity ratio on velxity profiles. 
Density ratio, 0.364; reference density, average of center- 
line density and coaxial-stream density. 




0 .4 .8 1.2 0 .4 .8 1.2 

Radial position, r 


(a) Velxity ratio, 10. (b) Velxity ratio, 30. 

Figure 7. - Effect of initial velxity ratio on mass fraction pro- 
files. Density ratio, 0.364; reference density, average of 
centerline density and coaxial-stream density. 
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Mass fraction of component 1, y Velocity, (u - u e )/(l - u 0 ) 




Radial position, r 

(a) Density ratio, 0. 364. <bl Density ratio, 0. 170. 

Figure 8. - Effect of initial density ratio on velocity profiles. 
Velxity ratio, 30; reference density, average of centerline 
density and coaxial -stream density. 




0 .4 .8 1.2 0 .4 .8 1.2 

Radial position, r 

(a) Density ratio, 0. 364. (b) Density ratio, 0. 170. 

Figure 9. - Effxt of initial density ratio on mass fraction profiles. 
Velxity ratio, 30; reference density, average of centerline den- 
sity and coaxial -stream density. 








0 .4 .8 1.2 0 .4 .8 1.2 

Radial position, r 


(a) Reference density, (b) Reference density, 

centerline density. coaxial-stream density. 

Figure 10. - Effect of reference density on velocity profiles. 
Velocity ratio, 20; density ratio, 0.182. 



0 .4 .8 1.2 0 .4 .8 1.2 

Radial position, r 


(a) Reference density, (b) Reference density, 

centerline density. coaxial-stream density. 

Figure 11. - Effect of reference density on mass fraction profiles. 
Velocity ratio, 20; density ratio, 0. 182. 
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Dimensionless mass of component 1, I 



0 1 2 3 4 5 


Axial position, x 

Figure 12. - Effect of initial velocity ratio on dimen- 
sionless mass of component 1. Density ratio, 
0.364; reference density, average of centerline 
density and coaxial -stream density; ratio of reac- 
tor radius to jet radius, 4. 



Figure 13. - Effect of initial density ratio on 
dimensionless mass of component 1. Ve- 
locity ratio, 30; reference density, aver- 
age of centerline density and coaxial- 
stream density; ratio of reactor radius to 
jet radius, 4. 
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